void Navier::CompHM1(int i, int j, int k, int gama)//calculate temporary W velocity
{
              //cout<<"\n ****************** Entrei na CompFV1 ******************************** \n";
              //cout<<"\n *****************************"<<i<<","<< j<<","<<k<<"***********************\n";
	      int alpha,beta;
	      double e0,e1,e2,c1,c2,d1,d2;
	      double coef0,coef1,coef2;
              double d2vdx2, d2vdy2, d2vdz2;
              double Gy;
              Gy = S.g[0]/S.froude;
              S.nu = 0.0;
	      alpha = 2*gama -1;
	      beta = alpha + gama;

	      coef0 = alpha*W[i][j][k+gama];
	      if(coef0>=0);
              {
		  coef0 *=(-1.0/DZ[j]);
		  e0= W[i][j][k+1-gama];
                  coef1 = U[i+1][j][k]+U[i][j][k];
		  coef2 = V[i][j+1][k]+V[i][j][k];
	      }
	     else
	     {
		  coef0 *=(1.0/DZ[j+alpha]);
	  	  if((k<S.gridp[2] || (alpha == -1)) 
		  {
		      e0 =.25*(gridVertexUVW.gridW(i,j,k+beta)+gridVertexUVW.gridW(i+1,j,k+beta) + 
			    gridVertexUVW.gridW(i,j+1,k+beta)+gridVertexUVW.gridW(i+1,j+1,k+beta));
		  }
		  else
		  {
		      e0 = W[i][j][k+1];/*Se o fluxo está vindo de uma ghost cell ela só pode ser inflow */
		  }

		  if(ISNOTAIR_OUTFLOW(flag[i-1][j][k+alpha]))
		  {
                     c1 = U[i][j][k+alpha];
		  }
                 else 
		 {
		     if(k<S.gridp[2] || (alpha == -1))
		    {	
                     	c1 = .25*(gridVertexUVW.gridU(i,j,k+gama)+gridVertexUVW.gridU(i,j,k+beta) + 
			   gridVertexUVW.gridU(i,j+1,k+gama)+gridVertexUVW.gridU(i,j+1,k+beta));
		    }
		    else
		    {
			c1 = .5*(gridVertexUVW.gridU(i,j,k+1)+ gridVertexUVW.gridU(i,j+1,k+1));
		    }
		 }
			
                 if(ISNOTAIR_OUTFLOW(flag[i+1][j][k+alpha]))
		 {
                     d1 = U[i+1][j][k+alpha];
		 }
                 else 
		 {
		     if(k<S.gridp[2] || (alpha == -1))
		     {
                     	d1 =.25*(gridVertexUVW.gridU(i+1,j,k+gama)+gridVertexUVW.gridU(i+1,j,k+beta) + 
			     gridVertexUVW.gridU(i+1,j+1,k+gama)+gridVertexUVW.gridU(i+1,j+1,k+beta));
		     }
		     else
		    {
			d1 = .5*(gridVertexUVW.gridV(i+1,j,k+1)+ gridVertexUVW.gridV(i+1,j+1,k+1));
		    }
		 }
		  coef1 = c1+d1;
		
		  if(ISNOTAIR_OUTFLOW(flag[i][j-1][k+alpha]))
		  {
                     c2 = V[i][j][k+alpha];
		  }
                 else 
		 {
		     if(k<S.gridp[2] || (alpha == -1))
		    {	
                     	c2 = .25*(gridVertexUVW.gridV(i,j,k+gama)+gridVertexUVW.gridV(i,j,k+beta) + 
			   gridVertexUVW.gridV(i+1,j,k+gama)+gridVertexUVW.gridV(i+1,j,k+beta));
		    }
		    else
		    {
			c2 = .5*(gridVertexUVW.gridV(i,j,k+1)+ gridVertexUVW.gridV(i+1,j,k+1));
		    }
		 }
			
                 if(ISNOTAIR_OUTFLOW(flag[i][j+1][k+alpha]))
		 {
                     d2 = V[i][j+1][k+alpha];
		 }
                 else 
		 {
		     if(k<S.gridp[2] || (alpha == -1))
		     {
                     	d2 =.25*(gridVertexUVW.gridV(i,j+1,k+gama)+gridVertexUVW.gridV(i,j+1,k+beta) + 
			     gridVertexUVW.gridV(i+1,j+1,k+gama)+gridVertexUVW.gridV(i+1,j+1,k+beta));
		     }
		     else
		    {
			d2 = .5*(gridVertexUVW.gridV(i,j+1,k+1)+ gridVertexUVW.gridV(i+1,j+1,k+1));
		    }
		 }
		  coef2 = c2 +d2;

	      }
	      if((coef1)>= 0)
	      {
		 coef1 *=(-1.0/(DX[i-1]+DX[i]));
                 if(ISFLUID_INFLOW(flag[i-1][j][k]))
                 {
                    f1=e1=W[i-1][j][k+gama];
		 }
		 else
		 {
		    f1=W[i][j][k+gama];
		    e1 =.25*(gridVertexUVW.gridW(i-1,j,k+gama)+gridVertexUVW.gridW(i-1,j+1,k+gama) + 
			     gridVertexUVW.gridW(i,j,k+gama)+gridVertexUVW.gridW(i,j+1,k+gama));
		 }
		 if(ISFLUID_INFLOW(flag[i+1][j][k]))
		 {
		    g1=W[i+1][j][k+gama];
		 }
		 else
		 {
		    g1=W[i][j][k+gama];
		 }
	       }
	       else
	       {
		  coef1 *=(1.0/(DX[i+1]+DX[i]))
                  if(ISFLUID_INFLOW(flag[i+1][j][k]))
                  {
                     g1=e1=W[i+1][j][k+gama];
		  }
		  else
		  {	
		     if(i<S.gridp[0])
		     {
			 g1=W[i][j][k+gama];
			 e1 =.25*(gridVertexUVW.gridW(i+1,j,k+gama)+gridVertexUVW.gridW(i+1,j+1,k+gama) + 
			     gridVertexUVW.gridW(i+2,j,k+gama)+gridVertexUVW.gridW(i+2,j,k+gama));
		     }
		     else
		     {
			 e1 =.5*(gridVertexUVW.gridW(i+1,j,k+1)+gridVertexUVW.gridW(i+1,j,k+1+1));
		     }
		   }
		   if(ISFLUID_INFLOW(flag[i-1][j][k]))
		   {
		      f1=W[i-1][j][k+gama];
		   }
		   else
		   {
		      f1=W[i][j][k+gama];
		   }
	      }
	      if((coef2)>=0)
	      {
		   coef2 *=(-1.0/(DY[j-1]+DY[j]))
                   if(ISFLUID_INFLOW(flag[i][j-1][k]))
                   {
                       f2=e2=W[i][j-1][k+gama];
		   }
		   else
		   {
		        f2=W[i][j][k+gama];
			e2 =.25*(gridVertexUVW.gridW(i,j-1,k+gama)+gridVertexUVW.gridW(i+1,j-1,k+gama) + 
			     gridVertexUVW.gridW(i,j,k+gama)+gridVertexUVW.gridW(i+1,j,k+gama));
		   }
		   if(ISFLUID_INFLOW(flag[i][j+1][k]))
		   {
			g2=W[i][j+1][k+gama];
		   }
		   else
		   {
		        g2=W[i][j][k+gama];
		   }
	       }
	       else
	       {
		  coef2 *=(1.0/(DY[j+1]+DY[j]))
                  if(ISFLUID_INFLOW(flag[i][j+1][k]))
                  {
                       g2=e2=W[i][j+1][k+gama];
		  }
		  else
		  {	
		       g2=W[i][j][k+gama];
		       if(j<S.gridp[1])
		       {
			   e2 =.25*(gridVertexUVW.gridW(i,j+1,k+gama)+gridVertexUVW.gridW(i+1,j+1,k+gama) + 
			        gridVertexUVW.gridW(i,j+2,k+gama)+gridVertexUVW.gridW(i+1,j+2,k+gama));
		       }
		       else
		       {
			   e2 =.5*(gridVertexUVW.gridW(i,j+1,k+1)+gridVertexUVW.gridW(i+1,j+1,k+1));
		       }
		  }
		   if(ISFLUID_INFLOW(flag[i][j+1][k]))
		   {
			f2=W[i][j-1][k+gama];
		   }
		   else
		   {
		        f2=W[i][j][k+gama];
		   }
	       }
	       double derconv= coef0 *(e0-W[i][j][k+gama])+ coef1*(e1-W[i][j][k+gama]) + coef2*(e2-W[i][j][k+gama]);

               d2wdz2=DD(W[i][j][k],W[i][j][k+gama],W[i][j][k+1],2,k+gama);
               d2wdx2=DDS(f1,W[i][j][k+gama],g1,0,i);
               d2wdy2=DDS(f2,W[i][j][k+gama],g2,1,j);

               H[i][j][k+gama] = S.nu*(d2wdx2+d2wdy2+d2wdz2) - derconv + Gz; 

               double temp=fabs(-derconv*S.delt +  W[i][j][k+gama]);
               double max_temp = 0;
               for(int idel=-1;idel<=1;idel++) 
               {
                  for(int jdel=-1;jdel<=1;jdel++)
                  {
                      for(int kdel=-1;kdel<=1;kdel++)
                      {
                          if((ISNOTAIR(flag[i+1+idel][j+jdel][k+gama+kdel])) &&
                             (fabs(W[i+idel][j+jdel][k+kgama+kdel]) > max_temp))
                                max_temp = fabs(W[i+idel][j+jdel][k+gama+kdel]);
                       }
                  }
               }
               if(temp > max_temp)
               {
                   cout<<"************************************************************************** \n";
                   cout<<"**********i+1,j,k,temp,max_temp"<<i<<","<<j<<","<<k<<","<<temp<<","<<max_temp<<"\n";
               }
}